classdef PlaneElement < Element
    % QUAD4ELEMENT 四节点等参单元
    % 用于平面应力、平面应变问题的四边形单元
    
    properties (Constant)
        ElementType = 'Quad4'
        NumNodes = 4
        Dimension = 2
        DofPerNode = 2
    end
    
    methods
        function obj = PlaneElement(connectivity, material, integrationRule)
            % 构造函数
            if nargin == 0
                connectivity = [];
            end

            if nargin <= 2
                integrationRule = GaussIntegrationRule(2); % 默认2×2高斯积分
            end
            
            obj@Element(connectivity, material, integrationRule);
        end
        
        function Ke = computeStiffnessMatrix(obj)
            % 计算单元刚度矩阵

            obj.validateConfiguration();
            
            nDofs = obj.getTotalDofs();
            Ke = zeros(nDofs, nDofs);
            
            % 获取积分点和权重
            [gaussPoints, weights] = obj.IntegrationRule.getPointsWeights();
            
            % 材料本构矩阵
            D = obj.MaterialModel.getConstitutiveMatrix();

            for i = 1:length(weights)
                xi = gaussPoints(i, 1);
                eta = gaussPoints(i, 2);
                weight = weights(i);
                
                % 计算应变-位移矩阵B和雅可比行列式
                natural_coord = [xi, eta];
                [B, detJ] = obj.computeStrainDisplacementMatrix(natural_coord);
                
                % 刚度矩阵贡献
                Ke = Ke + B' * D * B * detJ * weight;
            end
         
        end
        
        function Me = computeMassMatrix(obj)
            % 计算单元质量矩阵       
            obj.validateConfiguration();
            
            nDofs = obj.getTotalDofs();
            density = obj.MaterialModel.Density;
            
            if strcmpi(massType, 'lumped')
                Me = obj.computeLumpedMassMatrix(density);
            else
                Me = obj.computeConsistentMassMatrix(density);
            end
        end
        
        function Fe = computeLoadVector(obj, loadCase)
            % 计算单元载荷向量
            obj.validateConfiguration();
            
            nDofs = obj.getTotalDofs();
            Fe = zeros(nDofs, 1);
            
            switch lower(loadCase.Type)
                case 'bodyforce'
                    Fe = obj.computeBodyForceVector(loadCase.Value);
                case 'traction'
                    Fe = obj.computeTractionVector(loadCase);
                case 'pressure'
                    Fe = obj.computePressureVector(loadCase.Value);
                otherwise
                    error('Unsupported load type: %s', loadCase.Type);
            end
        end
        
        function stress = computeStress(obj, nodalDisplacements)
            % 计算单元应力
            obj.validateConfiguration();
            
            [gaussPoints, ~] = obj.IntegrationRule.getPointsWeights();
            nPoints = size(gaussPoints, 1);
            
            stress = struct();
            stress.gaussPoints = zeros(nPoints, 6); % [σxx, σyy, σzz, τxy, τyz, τzx]
            stress.vonMises = zeros(nPoints, 1);
            
            D = obj.MaterialModel.computeConstitutiveMatrix();
            
            for i = 1:nPoints
                xi = gaussPoints(i, 1);
                eta = gaussPoints(i, 2);
                
                % 计算应变-位移矩阵B
                [B, detJ] = obj.computeStrainDisplacementMatrix([xi, eta]);
                
                % 计算应变和应力
                strain = B * nodalDisplacements;
                stressPoint = D * strain;
                
                stress.gaussPoints(i, 1:4) = stressPoint;
                stress.vonMises(i) = obj.computeVonMisesStress(stressPoint);
            end
        end
        
        function strain = computeStrain(obj, nodalDisplacements)
            % 计算单元应变
            obj.validateConfiguration();
            
            [gaussPoints, ~] = obj.IntegrationRule.getPointsWeights();
            nPoints = size(gaussPoints, 1);
            
            strain = zeros(nPoints, 3); % [εxx, εyy, γxy]
            
            for i = 1:nPoints
                xi = gaussPoints(i, 1);
                eta = gaussPoints(i, 2);
                
                [B, detJ] = obj.computeStrainDisplacementMatrix([xi, eta]);
                strain(i, :) = B * nodalDisplacements;
            end
        end
        
        function N = evaluateShapeFunctions(obj, naturalCoords)
            % 计算形函数值
            xi = naturalCoords(1);
            eta = naturalCoords(2);
            
            N = zeros(1, 4);
            N(1) = 0.25 * (1 - xi) * (1 - eta);
            N(2) = 0.25 * (1 + xi) * (1 - eta);
            N(3) = 0.25 * (1 + xi) * (1 + eta);
            N(4) = 0.25 * (1 - xi) * (1 + eta);
        end
        
        function dN_dxi = evaluateShapeGradients(obj, naturalCoords)
            % 计算形函数导数
            xi = naturalCoords(1);
            eta = naturalCoords(2);
            
            % 自然坐标下的形函数导数
            dN_dxi = zeros(2, 4);
            
            dN_dxi(1, 1) = -0.25 * (1 - eta);
            dN_dxi(1, 2) =  0.25 * (1 - eta);
            dN_dxi(1, 3) =  0.25 * (1 + eta);
            dN_dxi(1, 4) = -0.25 * (1 + eta);
            
            dN_dxi(2, 1) = -0.25 * (1 - xi);
            dN_dxi(2, 2) = -0.25 * (1 + xi);
            dN_dxi(2, 3) =  0.25 * (1 + xi);
            dN_dxi(2, 4) =  0.25 * (1 - xi);
        end
        
        function J = computeJacobian(obj, dN_dxi)
            % 计算雅可比矩阵 
            J = zeros(2, 2);
            for i = 1:4
                J(1,1) = J(1,1) + dN_dxi(1,i) * obj.NodeCoordinates(i, 1);
                J(1,2) = J(1,2) + dN_dxi(1,i) * obj.NodeCoordinates(i, 2);
                J(2,1) = J(2,1) + dN_dxi(2,i) * obj.NodeCoordinates(i, 1);
                J(2,2) = J(2,2) + dN_dxi(2,i) * obj.NodeCoordinates(i, 2);
            end
        end
     
    end
    
    methods (Access = protected)
        function [B, detJ] = assembleStrainDisplacementMatrix(obj, naturalCoords)
            % 组装应变-位移矩阵B
            dN_dxi = obj.evaluateShapeGradients(naturalCoords);
         
            J = obj.computeJacobian(dN_dxi); % 使用缓存的dN_dxi
            detJ = det(J);
            % 计算物理坐标下的形函数导数
            dN_dx = J \ dN_dxi;
            
            B = zeros(3, 8); % 平面问题：3个应变分量，8个自由度
            
            for i = 1:4
                % B矩阵组装
                B(1, 2*i-1) = dN_dx(1, i); % εxx
                B(2, 2*i)   = dN_dx(2, i); % εyy
                B(3, 2*i-1) = dN_dx(2, i); % γxy
                B(3, 2*i)   = dN_dx(1, i);
            end
        end
        
        function Me = computeConsistentMassMatrix(obj, density)
            % 计算一致质量矩阵
            nDofs = obj.getTotalDofs();
            Me = zeros(nDofs, nDofs);
            
            [gaussPoints, weights] = obj.IntegrationRule.getPointsWeights();
            
            for i = 1:length(weights)
                xi = gaussPoints(i, 1);
                eta = gaussPoints(i, 2);
                weight = weights(i);
                
                N = obj.evaluateShapeFunctions([xi, eta]);
                [~, detJ] = obj.computeJacobian([xi, eta]);
                
                % 组装N矩阵 (2×8)
                N_mat = zeros(2, 8);
                for j = 1:4
                    N_mat(1, 2*j-1) = N(j);
                    N_mat(2, 2*j)   = N(j);
                end
                
                Me = Me + N_mat' * N_mat * density * detJ * weight;
            end
        end
        
        function Me = computeLumpedMassMatrix(obj, density)
            % 计算集中质量矩阵
            nDofs = obj.getTotalDofs();
            Me = zeros(nDofs, nDofs);
            
            volume = obj.computeVolume();
            totalMass = density * volume;
            nodalMass = totalMass / obj.NumNodes;
            
            % 对角线元素赋值
            for i = 1:obj.NumNodes
                dofIndices = (2*i-1):(2*i);
                Me(dofIndices, dofIndices) = nodalMass * eye(2);
            end
        end
        
        function Fe = computeBodyForceVector(obj, bodyForce)
            % 计算体力载荷向量
            nDofs = obj.getTotalDofs();
            Fe = zeros(nDofs, 1);
            
            [gaussPoints, weights] = obj.IntegrationRule.getPointsWeights();
            
            for i = 1:length(weights)
                xi = gaussPoints(i, 1);
                eta = gaussPoints(i, 2);
                weight = weights(i);
                
                N = obj.evaluateShapeFunctions([xi, eta]);
                [~, detJ] = obj.computeJacobian([xi, eta]);
                
                % 组装N矩阵
                N_mat = zeros(2, 8);
                for j = 1:4
                    N_mat(1, 2*j-1) = N(j);
                    N_mat(2, 2*j)   = N(j);
                end
                
                Fe = Fe + N_mat' * bodyForce * detJ * weight;
            end
        end
        
        function Fe = computeTractionVector(obj, tractionCase)
            % 计算面力载荷向量
            % 简化实现，假设面力作用在单元边上
            Fe = zeros(obj.getTotalDofs(), 1);
            % 具体实现需要根据面力作用的边来积分
        end
        
        function Fe = computePressureVector(obj, pressure)
            % 计算压力载荷向量
            Fe = zeros(obj.getTotalDofs(), 1);
            % 具体实现需要根据压力方向和作用面来积分
        end
        
        function vonMises = computeVonMisesStress(obj, stress)
            % 计算von Mises等效应力
            if length(stress) == 3 % 平面应力
                sxx = stress(1);
                syy = stress(2);
                txy = stress(3);
                vonMises = sqrt(sxx^2 - sxx*syy + syy^2 + 3*txy^2);
            else
                vonMises = 0;
            end
        end
        
        function validateConfiguration(obj)
            % 验证单元配置
            if isempty(obj.NodeCoordinates)
                error('Node coordinates must be set before computation');
            end
            if isempty(obj.MaterialModel)
                error('Material model must be set before computation');
            end
            if isempty(obj.IntegrationRule)
                error('Integration rule must be set before computation');
            end
        end
    end
end